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Abstract 



In this thesis we present a kinetic Monte Carlo model for the description of epitaxial graphene 
growth. Experimental results suggest a growth mechanism by which clusters of 5 carbon atoms 
are an intermediate species necessary for nucleation and island growth. This model is proposed 
by experimentally studying the velocity of growth of islands which is a highly nonlinear function 
of adatom concentration. In our simulation we incorporate this intermediate species and show 
that it can explain all other experimental observations: the temperature dependence of the 
adatom nucleation density, the equilibrium adatom density and the temperature dependence of 
the equilibrium island density. All these processes are described only by the kinematics of the 
system. 
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Chapter 1 

Introduction 



1.1 Scope of this thesis 

Graphene is a 2-dimensional layer of graphite, formed by carbon atoms arranged in an hexagonal 
lattice. For a long time, graphene was used as the starting point for the theoretical study of 
carbon allotropes, such as carbon nanotubes, which can be described by the rolling of a graphene 
sheet into the form of a cylinder. But it was thought that graphene was merely a theoretical 
construct because thermal fluctuations would render the isolated 2-dimensional sheets unstable. 
Nonetheless, it was reported by Novoselov and co-workers in 2004 [lj that graphene sheets had 
been isolated and some of their physical properties measured. Since then, it has attracted the 
attention of the physics and materials communities [2] due to its interesting physical properties 
and its potential for technological applications. The 2010 Physics Nobel Prize was awarded to 
Andre Geim and Konstantin Novoselov 

for groundbreaking experiments regarding the two-dimensional material graphene 

as described by the awarding body. 

Recent efforts have focused on synthesising graphene sheets by means of epitaxial growth, 
one of the most promising routes to large scale and high quality graphene production. There 
exist ample studies characterising the properties of graphene sheets grown on a variety of sub- 
strates. However, very little work has been done on the growth kinetics of the graphene sheets. 
The first experimental studies only became available in 2008 [31 d] . The remarkable behaviour 
observed in these studies is very different from that of known metal-on-metal growth systems. 
Initial theoretical studies inspired by the experimental findings have been conducted using a 
rate equations approach [5], and these studies have been able to describe the main features of 
the experimental results. 
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In this thesis we present a kinetic Monte Carlo model for theoretical studies of epitaxial 
graphene growth kinetics. The model is based on the basic principles found experimentally 
and used in the rate equations study, and provides further insight into the physical processes 
governing the system. 

1.2 Thesis outline 

Chapter |2j is organised in two parts, and it describes some background information relevant to 
the project. First we describe the physical properties and technological potential of graphene, 
which have made it an attractive system for both the physics and materials communities. In 
the second part we present epitaxial growth, discussing the physical processes involved together 
with the experimental techniques used in the field. 

In Chapter [3] we give a detailed overview of the current understanding of epitaxial graphene 
growth. We survey the most relevant experimental and theoretical results describing the kinetics 
of epitaxial graphene growth, and abstract the principal conclusions to be incorporated in a 
model to study the system further. 

In Chapter [4] we introduce in detail the kinetic Monte Carlo (KMC) technique, focusing on 
the N-fold way algorithm. We also provide motivation for the use of this method in the study 
of epitaxial growth systems. 

In Chapter [5] we present a simple standard KMC model for epitaxial growth, exemplifying 
the N-fold way algorithm. This simulation will be used as a starting point for more complex 
models appropriate for the study of the graphene system. 

In Chapter [6] we present a KMC model that includes the main features observed in graphene 
growth experiments. All the experimental observations are explained physically within this 
model, identifying the most important processes governing the kinetics of epitaxial graphene 
growth. 

In Chapter [7] we summarise the most relevant results and discuss lines of research that can 
be followed from the present work. 



2 



Chapter 2 

Graphene and epitaxial growth 



In this chapter we first present the most relevant physical properties of graphene and discuss 
its technological potential. We then look at both theoretical foundations and experimental 
techniques in epitaxial growth, which is most probably the method of choice for large scale and 
high quality graphene production. 



2.1 Graphene 

2.1.1 Physical properties 

The physical properties of graphene were recently reviewed in Ref. [6] . Graphene is a 2-dimensional 
layer formed by sp 2 -hybridised carbon atoms arranged in a hexagonal lattice and separated by 
a ~ 1.42 A. The hybrid orbitals give rise to a bonds between atoms. The remaining orbital, 
p z , is perpendicular to the planar structure and forms covalent bonds between neighbouring 
carbon atoms, leading to a half filled tt band. This band can be described using a tight-binding 
approach with a single hopping matrix element between neighbouring atoms —t. The resulting 
band structure, first calculated by Wallace in 1947 [7], has a dispersion relation 



E(k) = ±t^ 3 + 2 cos(V3k y a) + 4 cos ^-^-/^a^ cos (^k x a^j (2.1) 
for wavevector (k x ,k y ). 

Graphene is a zero gap semiconductor, where the Fermi surface consists only of six Fermi 



points at the edge of the Brillouin zone, as shown in Fig. 2.1 Expanding about the Fermi points 
kp = k — q in reciprocal space we find E±(q) ~ drup|q| + 0[(q/kp) 2 ] for constant vp ~ 10 6 
ms -1 . This parallels the dispersion relation of ultrarelativistic particles described by the massless 
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Figure 2.1: Electronic dispersion relation of graphene obtained with a tight-binding calculation with nearest 
neighbour interactions, given by Eq.(2.1l. The energy -E(k) is in units of t. 



Dirac equation, so the Fermi points are commonly called Dirac points. It is this dispersion that 
determines most of the singular physical properties of graphene and it means that graphene is 
a laboratory condensed matter system to test (2+l)-dimensional quantum electrodynamics [8]. 

As reported in Ref . [U] the experimental study of some of the physical properties confirmed the 
existence of massless Dirac carriers in graphene. For instance, their cyclotron mass depends on 
the square root of the density of states, and the integer quantum Hall effect occurs at half-integer 
filling factors, both characteristic of massless Dirac fermion systems. 

2.1.2 Technological applications 

The physical properties of graphene result in large carrier mobilities that persist at room tem- 
perature, even with the presence of doping species. This means graphene holds the potential to 
become a replacement for silicon in the electronics industry [2J, along with a wider variety of 
applications. 

Transistors and diodes need the presence of a band-gap for their operation, so standard 
graphene sheets are not appropriate. Researchers have explored alternative graphene-based 
structures with the presence of a band-gap, for instance graphene nanoribbons |10|lll] where the 
band-gap is proportional to the ribbon width, or graphene nanomeshes [12] where the graphene 
sheets are punched with an array of nanoscale holes. Graphene-made transistors in the GHz 
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scale were recently reported by IBM researchers \1'6\ [T^l [To] with performances superior to those 
of similar silicon transistors. 

Other examples are studies on the large heat conductivity of graphene [HI [T7] with potential 
applications in nanoelectronics where large heat dissipation is needed. 

2.2 Epitaxial growth 

Epitaxial growth is the name given to the process of producing epitaxial thin films on substrates. 
It is a widespread technique with applications ranging from the production of semiconducting 
devices to nanotechnology and it has also attracted the scientific community because of the 
complex atomic processes involved. 

2.2.1 Theoretical description 

Epitaxial growth can be classified in three so-called modes, first introduced in the seminal 
work by Bauer |18j . Following Ref.[19], the different growth modes can be understood by 
thermodynamic arguments. In the layer, or Frank-van der Merwe mode, the atoms are more 
strongly attracted to the substrate than to themselves, and the epitaxial film is formed layer 
after layer. Quantitatively, it arises in the deposition of material A on B when the surface free 
energies ji for i = A,B obey 

1A + lint < IB (2.2) 

for interface free energy 7^ between surfaces A and B. This follows because the free energy 
needs to be minimised at equilibrium, and the inequality requires to maximise the area covered 
by deposit A. The island, or Volmer-Weber mode, results when the atoms are more strongly 
attracted by each other than to the substrate, and multilayer (3-dimensional) islands form. 
Quantitatively, 7^4 + ji n t > 1b- A third hybrid mode termed Stranski-Krastanov growth mode, 
or layer-plus-island, arises because the interface energy 7^ increases as the layer thickness 
increases, so island growth starts as a layer mode but turns into an island mode. 

In the present work we only consider submonolayer growth when the first layer is forming, 
because we only reach fractional coverages of the lattice and the system is purely 2-dimensional. 

In the atomic regime there are many processes during deposition and growth. Atoms are 
deposited at a certain rate on the substrate, and then undergo a series of processes. They 
can re-evaporate, diffuse over the substrate or along island edges, nucleate to form islands, join 
growing islands or detach from islands. Temperature is a key parameter because the different 
atomic processes are thermally activated. In the present work we incorporate these atomistic 
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processes into a kinetic Monte Carlo model. 
2.2.2 Experimental techniques 

Molecular beam epitaxy (MBE) [2D] is the most widespread technique used in epitaxial growth. 
A beam of atoms or molecules is deposited on a previously prepared substrate kept at high 
temperatures to allow the arriving particles to diffuse over its surface. The deposition is carried 
out under ultra-high vacuum conditions in order to minimise impurities. Control over the beam 
allows films grown using MBE to be of very high quality and to have the desired properties. 
Typical deposition rates are ~1 ML/s, which are high enough to significantly reduce the incor- 
poration of impurities into the growing material. Another technique used in epitaxial growth is 
vapour phase epitaxy where the substrate is placed in contact with a gas containing the deposit 
elements, and reactions between the two lead to epitaxial growth. 

The grown sheets can be observed with high accuracy using different microscopy techniques 
such as scanning tunneling microscopy or atomic force microscopy. Furthermore, the growing 
process can be monitored using low energy electron microscopy (LEEM). For example, in some 
experiments [31 0] relevant for the work described in this thesis the changes in reflectivity of a 
LEEM are used to infer the carbon adatom concentration on the substrate. 
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Chapter 3 

Experimental and theoretical results 
on epitaxial graphene growth 

Epitaxial growth is one of the most promising techniques currently being explored to synthe- 
sise high quality graphene sheets with the appropriate properties for technological applications. 
Much ongoing research concentrates in the production of such epitaxial sheets and their char- 
acterisation for a variety of substrates with different physical properties have been reported. 
However, studies of the growth kinetics of graphene are just beginning both experimentally and 
theoretically and here we present a review of the current status of such studies. 



3.1 Experimental results 

In this section we present the most relevant experimental results that motivated the present 
work. They are the first that look at the kinetics of the growth of epitaxial graphene. 



3.1.1 Kinetics of epitaxial graphene growth 

The most relevant experimental observations of the kinetics of epitaxial graphene growth have 
been reported by Loginova and co-workers in Refs.[3l 0]. They deposit pure carbon atoms 
and ethylene molecules on ruthenium Ru(0001) and iridium Ir(lll) substrates under ultra-high 
vaccum conditions and high substrate temperatures. Low Energy Electron Microscopy (LEEM) 
is used to monitor the growth of graphene islands and the local carbon adatom concentration 
on the substrate. 

A typical adatom concentration profile that will appear throughout this work is as shown 



in Fig. 3.1 A deposition flux F is turned on at about time t = 150 s. Initially the adatom 
concentration increases almost linearly with time due to the constant flux, in the example given 
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Figure 3.1: Typical adatom density n curve as a function of time found in epitaxial growth experiments, adapted 
from Loginova et al. [3]. 



of F = 0.0023 ML/s, small for typical epitaxial experiments. At about t = 250 s the adatom 
density curve reaches a maximum, called n nuc because it approximately corresponds to the 
onset of nucleation. Islands start nucleating and adatoms start disappearing from the substrate 
because they attach to islands. This decrease continues for some 100 s until an equilibrium 
between islands and adatoms is reached. Later on, at t = 700 s, the deposition flux is turned off, 
quickly leading to a different equilibrium concentration between adatoms and islands, smaller 
than the previous equilibrium because adatoms are no longer being incorporated externally into 
the system. The adatom density in this regime is labelled n e q. 

Note the definition of n nuc can be slightly confusing in systems where the critical island size 
is not well-defined. The critical island size i is an island size above which adatom attachment 
can be considered as irreversible, so that nucleation can be defined as the process of an island 
going from size i —> (i + 1). Then, in systems with no critical island size, large clusters that 
eventually lead to islands can start appearing before the maximum of the adatom density profile 
is reached. Throughout this work we are going to define n nuc as the maximum of the adatom 
density profile. 

In their studies, Loginova and co-workers find that the island growth velocity presents a 



nonlinear dependence on the carbon adatom concentration as shown in Fig. 3.2 The island 
growth velocity v is defined as v = P~ 1 dA/dt for island perimeter P and area A. In most 
growth systems the island growth velocity is found to be proportional to the supersaturation of 
adatoms, i.e. proportional to the difference in adatom concentration n and adatom equilibrium 



concentration n eq , v 



n, 



eq, 



for constant C. To explain the nonlinear relationship found 



for graphene growth, Loginova and co-workers propose a model in which the energy barrier for 
monomer attachment to islands is larger than the barrier for the formation of clusters of m 
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Figure 3.2: Experimental observation of nonlinear island growth velocity as a function of adatom concentration. 
In the plot, c nucl corresponds to n nuc in the text. Taken from Loginova et al. [3]. 



carbon atoms and their posterior attachment to islands. They assume that the growth velocity 
is proportional to the supersaturation of clusters rather than adatoms. Following Ref.[3], the 
concentration of m-clusters in a supersaturated adatom sea has an exponential dependence 
on the energy difference between m isolated carbon atoms and the energy needed to form an 
m— atom cluster E„ 



■"mi 



c (m) = e ( mfJ ,-E m )/k B T = I n \ e -E m /k B T ^ (g-q 



where the sea of carbon adatoms is assumed to be an ideal lattice gas with carbon chemical 
potential fj, = kBTln(n/n eq ). The island growth velocity as a function of adatom concentration 
is then 



OJc (m) -cth = B 



eq 



n , 

1 



(3.2) 



where B = C m e~ Em ^ kBT for C m the proportionality constant in the velocity dependence in 
cluster supersaturation. Fitting the data in Fig. |3.2[ Loginova and co-workers find a best 
estimate of m ~ 5. This means that in their model clusters formed by 5 carbon atoms are an 
intermediate species that determines the growth kinetics of epitaxial graphene growth. This 
intermediate species could be formed on the graphene free lattice and diffuse over the substrate 
to attach to islands as is assumed throughout Ref.[3|. This model is the one we are going to take 
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in our KMC approach to the description of the system. However, it is pointed out by Loginova 
and co-workers that the 5-clusters could instead form only near island edges at the moment of 
attachment to graphene. 
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Figure 3.3: Adatom density at nucleation n nuc ( c ' ItlcI i n the plot) and at equilibrium n eq (c eq in the plot) as a 
function of temperature. The data corresponds to deposition of carbon (filled labels) or ethylene (emptly labels). 
Taken from Loginova et al. [3]. 




Further observations in the graphene growth system that differ from other known growth 
scenarios are the temperature dependence of adatom density at nucleation and the adatom 
density at equilibrium. 



It can be seen in Fig. 3.3 that the nucleation density increases with increasing temperature. 
Most growth systems are diffusion limited, so that increasing the temperature results in higher 
adatom mobility and earlier nucleation. Therefore, the nucleation density usually decreases with 
increasing temperature unlike for graphene. 



It can also be observed in Fig. 3.3 that the nucleation concentration is roughly twice as 
large as the equilibrium density n nuc ~ 2n eq . This indicates a small energy barrier to adatom 
detachment from graphene sheets and a correspondingly large attachment barrier for adatoms 
to growing graphene. This observation is in agreement with the model by which the dominant 
species in the kinetics of graphene growth is an intermediate 5-cluster. 

Another experimental observation that remains unpublished is the behaviour of the island 
density at equilibrium. In agreement with known growth systems (described below in Chapter [5]) 
the island density decreases with increasing temperature. However, the decrease in the graphene 
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system is found to be much larger than in other systems, suggesting the existence of a mechanism 
that allows the occurence of a large number of nucleations at low temperatures. We thank Elena 
(Loginova) Starodub for providing the data related to island density. 

We finally note that the growth mechanism of graphene is different for other metal sub- 
strates. For instance, for Ru(0001) and Ir(lll), the deposition of molecules containing hydrogen 
and carbon decompose rapidly and growth is determined by individual carbon atoms. How- 
ever, theoretical and experimental studies |21| [22] indicate that graphene growth on copper is 
determined by the deposited molecules, and that only at late stages of the graphene formation 
process is hydrogen released from the composites. 



3.1.2 Graphene nanoclusters 

Another set of experiments |23[ [M] looks more carefully at the very initial steps of graphene 
growth when nucleation occurs. The studies are both on Ru(0001) and Ir(lll) as the experiments 
by Loginova and co-workers reported above. Small carbon nanoclusters of sizes of the order 
of tens of carbon atoms are observed on both substrates. They are described as a possible 
predecessor for graphene islands, and could be understood as a critical island size for graphene. 
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(c) Coronene (£) A-7C6 (l) B-3C6 

Figure 3.4: Images of coronene and dome-like carbon nanoclusters on Ru(0001), together with a depiction of 
the structure of the clusters. Taken from Cui et al. [23] . 



In both cases, the observed nanoclusters have a dome-like shape where the adatoms at the 
perimeter of the quasi-circular islands are strongly attached to the substrate and the adatoms 
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in the center are highly detached from it. This structure can be seen in Fig. 3.4 



Other studies of epitaxial graphene growth on Rhodium Rh(lll) reported in Ref.[25j also 
indicate the presence of these nanoclusters. The study indicates that the nanoclusters are mobile 
on Rh(lll), so they would not correspond to immobile critical islands. However, the experiment 
is performed on a different substrate to the above experiments, and the experimental method 
differs as well. This means that the results might not be directly comparable to the above. 



3.2 Theoretical results 

Based on the experiments reported by Loginova and co-workers in Refs.0 H], Zangwill and 
Vvedensky proposed a rate equations (RE) model [5] to describe epitaxial graphene growth. RE 
are described in the next chapter. They incorporate the intermediate 5-clusters and propose a 
system evolution determined by a set of coupled differential equations for the adatom density 
n, the 5-atom cluster density c and the island density N as follows 



dvt 

— = F — iDn 1 + iKc — DnN + K'N, (3.3) 
dc 

— = Dn l -Kc- D'cN - jD'c 3 , (3.4) 

dN . , . . 

- = DW, (3.5) 

The parameters are carbon atom deposition flux F, adatom diffusion rate D, cluster diffusion 
rate D', cluster dissolution rate K and adatom detachment rate from islands K' , all assumed to 
be of the Arrhenius form because they are thermally activated. The index i = 5 is the cluster 
size, and the index j is an unknown for the number of clusters that need to come together to 
form an island. 

The above set of equations together with an appropriate choice of values for the parameters 



reproduces the main features of the adatom density curve in Ref.jl] as shown in Fig. 3.5. The 
temperature dependence of the adatom density at nucleation is well-reproduced within the RE 
by introducing the variable j and taking a value j > 6, and with the presence of the cluster 
dissociation rate. However, the temperature dependence of island density at equilibrium cannot 
be explained with the above RE model. Discrepancies with the experimental data are attributed 
to spatial effects that are not included in the RE approach. 
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Figure 3.5: Rate equations comparison with LEEM data from Ref. [3] for the adatom density profile. Adapted 
from Zangwill and Vvedensky pj]. 

3.3 Kinetic Monte Carlo model 

Based on the above experimental and theoretical results of graphene growth, a KMC model to 
describe epitaxial growth of graphene sheets should incorporate the following basic ingredients: 
three different species O Hj corresponding to the experimentally observed carbon monomers, 
carbon 5-clusters, and graphene sheets; and a large critical nucleus size [23} [21]. 

These basic components of the model should form the basis for the description of the exper- 
imental observations on the kinetics of epitaxial graphene growth [31 0] : 

1. The adatom density at the onset of nucleation increases with increasing temperature. 

2. The adatom density at equibrium n eq is roughly half of the adatom density at the onset 
of nucleation n nuc , namely, n nuc ~ 2n eq . 

3. The island density at equilibrium decreases rapidly with increasing temperature. 

In the rest of this thesis we will present such a KMC model and describe the relevant physical 
processes in epitaxial graphene growth within the model. 
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Chapter 4 



Kinetic Monte Carlo 



In this chapter we present the kinetic Monte Carlo (KMC) method as used in the study of 
epitaxial growth. The most widespread algorithm used in the field is the so-called N-fold way 
algorithm, which is adopted in this thesis. 



4.1 Theoretical studies and computer simulations 

In the field of many-body dynamics simulations, and in particular epitaxial growth systems, 
there are various approaches with different levels of approximation and capabilities. 

A completely deterministic mean-field approach called rate equations (RE) is commonly used 
in studies of epitaxial growth. RE are a finite set of coupled first order differential equations 
for the densities of the different species in an epitaxial system and with a rate associated with 
each possible process. The strength of the technique is its simplicity, but the pay-off is in 
accuracy. Standard RE cannot include any spatial information, but more complex approaches 
using capture numbers can encapsulate it. Zangwill and Vvedensky have used RE to describe 



the epitaxial graphene growth system |5J as described above in Section 3.2 

Molecular dynamics (MD) simulations are used for detailed analysis of the dynamics of a 
system. The forces between different components are calculated and the system is evolved 
according to these. Both classical and quantum approaches to MD simulations are used, and 
major advances in the field, for instance the Car-Parrinello method [26J incorporating density 
functional theory into MD, have led to widespread use of the method. However, the timescales 
and system sizes that can be explored remain small, restricting the applicability of the technique. 

KMC looks at intermediate time and length scales and has proved appropriate for the de- 
scription of epitaxial growth systems. 
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4.2 Kinetic Monte Carlo 



In this section we describe the KMC technique. We first present the N-fold way algorithm used 
in this thesis and next we discuss why the KMC method can be used in the description of 
epitaxial growth processes. 



4.2.1 N-fold way algorithm 

The KMC method is based on the use of random numbers to simulate the time evolution of a 
system in discrete steps, in such a way that each step represents a move of the entire system from 
one state to another. We are going to concentrate on the N-fold way algorithm, first introduced 
by Bortz, Kalos and Lebowitz in 1975 |27j . 

Following Ref. [28J, a rate ri is associated with process i, and rates are usually taken to be 
of the Arrhenius form because they are thermally activated, 

Ti = v e- Ei / kBT (4.1) 

where vq = 2ksT/h is of the order of the atomic vibrational frequency, and Ei is the energy 
barrier associated with process i. T is the temperature, ks is Boltzmann's constant and h is 
Planck's constant. This form for the rates will imply in a KMC model that events that are more 
likely to occur will happen more frequently. At each time step, all rates are calculated for the 
given configuration of the system. 

To select the process that is going to occur at a given time step, a total transition rate R is 
constructed as 

JV 

fl = J> ( 4 - 2 ) 
i=i 

for a total of N possible processes in a given configuration of the system. A number p\ £ [0, 1) 
uniformly distributed is calculated, and then the process j such that 



j-l 3 

i=i i=i 



r l <Rpi<Y^r l (4.3) 



is selected to take place in the given time step. This is shown in Fig. 4.1 

Practically, possible events can be grouped together and the search time for the event is 
reduced. For instance, all free adatoms will have equal probability of diffusing into any of the 
four nearest-neighbour sites in a square lattice simulation, and all these events can be put into a 
subgroup. Constructing these subgroups reduces the total rate to a set of partial sums, and the 
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Rpi for random p\ £ [0, 1) 



Figure 4.1: Schematic of the N-fold way algorithm selecting the rate Ti to occur in the given time step. The 
size of the boxes containing the rates is proportional to the rates themselves. 



search then proceeds first by identifying the relevant subgroup and then searching only within 
this subgroup. This search method is referred to as the binning method [28] or 2-level scheme as 
the search is done in two levels. There are methods available |29j which are up to 7 times faster 
than the binning method and are based on a further subdivision of the list of possible events, 
refered to as if-level schemes for K subdivisions. 

The time evolution of the system is based on the assumption that the probability of an event 
occuring is independent of the history of the system and therefore it obeys Poisson statistics. 



This assumption will be justified below in Subsection 4.2.2| Following Ref.|27j. the probability 
that an event i occurs in the infinitessimal time interval (t, t + dt) is pi(t)dt = ridt where ri 
is the rate associated with the event. The total probability P(t) that an event occurs in time 
interval (t, t + dt) is then 

P(t)dt = Pi(t)dt = ^2 ndt = Rdt. (4.4) 

i i 

Let ir(t) be the probability that no event occurs in time interval (0,t). Then, the probability 
that no event occurs in time interval (0, t + dt) is equal to the probability that no event has 
occured in time interval (0, t) and the probability than no event occurs in time interval (t,t + dt), 
namely, 

ir(t + dt) = vr(t)(l - Rdt) = ir(t) - ir(t)Rdt. (4.5) 
After rearranging and taking the limit dt — > 0, we obtain 

= Hm *(*+<*;)-*(*) = _ <t)R (46) 

dt dt^o dt w v ; 
which can be easily solved to give 

?r(t) = Tr(0)e- Rt (4.7) 

with 7r(0) = 1 as the probability of no event occuring in time t = is unity. The physical time 



of each time step in the simulation is then assumed to be distributed as Eq.(4.7) and a second 
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random number p2 G (0,1) is chosen such that the physical time r is given via P2 = e Rt or 



This expression gives an average time between events equal to 1/R. 

4.2.2 Use of kinetic Monte Carlo models to describe epitaxial growth 

In this section we motivate the use of a stochastic Monte Carlo technique for the description of 
epitaxial growth systems. 

As described in Ref . |30j . deposition often occurs by the random impingement of atoms or 
molecules on a substrate. Within a KMC simulation, a deposition flux is included in the model, 
and the deposition sites are chosen randomly, in accordance with experimental realisations. 

After deposition, atoms are adsorbed on specific sites on the lattice that correspond to 
minima in energy for the configuration of the system. A given configuration of the adatoms on 
the substrate is called a state of the system. The adatoms vibrate at frequencies of the order 
uq ~ 10 13 s _1 about the minima they occupy. However, these vibrations do not change the 
state of the system, only the state of the individual adatoms. The state of the system changes 
occasionally when an adatom vibration makes the adatom move from one minimum to another 
by going past an energy barrier E. Such transitions are thermally activated, so the rate r 
associated with them is of the Arrhenius form. 

A KMC model to describe epitaxial growth takes advantage of the different time scales 
between the change of state of an adatom and of the entire system by only considering the 
latter. During the irrelevant vibrations of an adatom about a given minimum, the system 
loses memory about the previous configurations it occupied, and therefore one can use Poisson 
statistics to describe a KMC time step. 



r = — 



HP2) 
R 




18 



Chapter 5 



Standard kinetic Monte Carlo model 
of epitaxial growth 

In this chapter we introduce what we call the KMC standard model of epitaxial growth. The 
objective is two-fold: first, to provide an example of the N-fold way algorithm described above; 
and second, to have a simple model that is going to serve as a starting point for further studies 
and as a reference to compare with the graphene system. 

5.1 Description 

We consider a 200x200 square lattice with periodic boundary conditions. The model describes 
two processes: deposition of single atoms and diffusion of adatoms (adsorbed atoms) over the 
substrate. This model does not aim to describe the graphene system, but the studies we will 
perform will be keeping in mind the final objective to study epitaxial graphene growth. 

Deposition occurs at a constant flux F, and atoms can only be deposited at an unoccupied 
lattice location. In the case that a deposition attempt is made on an occupied site, then a 
deposition process is still forced to happen in order to keep the deposition flux constant, and a 
different site is chosen repeatedly until an unoccupied site is found. The total coverages explored 
here are a small fraction of the lattice because we are interested in the initial stages of epitaxial 
graphene growth. Then, such events are rare indicating that the above modelling of deposition 
does not introduce an undesired effect in the model. 

Diffusion occurs between nearest-neighbour sites, and an adatom can only diffuse to an 
unoccupied site. The diffusion rate D is given by the Arrhenius form 

D = Uoe -(E a +nE b )/k B T (5.1) 
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where vq = 2kgT /h for Boltzmann constant fee, temperature T and Planck constant h. Fur- 
thermore, E a is the energy barrier associated with the diffusion of an adatom, and E\, is an extra 
energy barrier associated with a bond to a nearest neighbour. The variable n is the number 
of nearest neighbours of the diffusing adatom, and the above bond-counting model is typically 
used in epitaxial growth simulations [30J. Note this rate depends on the local configuration 
around the diffusing adatom, where nearest-neighbours result in a penalty energy for diffusion. 
That is, adatoms that share bonds with others are less likely to diffuse away and this is the key 
ingredient to get island formation. 

Temperature is the parameter we tune in order to explore the different regimes of the system. 
The energy barrier for adatom diffusion is kept fixed at E a = 1.00 eV, while the bond energy 
barrier Eb is varied. The flux is set to F = 1 ML/s and the total coverage reached is 9 = 0.1 
ML. The units of coverage are monolayers (ML) where one monolayer corresponds to a layer of 
epitaxial material grown on the substrate. 

The model presented here belongs to a class of models referred to as generic models |30j . 
These generic models do not have any specifications of critical island size, that is, no island 
size above which attachment is irreversible is defined. Instead, a set of local rules such as the 
bond-counting scheme above is defined and the system is left to evolve according to these only. 



5.2 Results and discussion 

In this section we analyse the physical processes involved in epitaxial growth based on the KMC 
model described above. We start with a qualitative look at the growth behaviour, and then we 
move on to a quantitative analysis that will serve as a basis for the later study of the graphene 
system. 



5.2.1 Physical lattice 



In Fig. |5.1| we can see the physical lattice for different temperatures. All instances are after 
the last adatom of a total coverage of 9 = 0.1 ML has been deposited, corresponding to a 10% 
coverage of the lattice. The bond energy barrier used is E^ = 0.20 eV. 



It can be seen from Fig. 5.1 that as temperature increases the islands that are formed 
become larger in size and fewer in number. Adatoms at the edges of islands are more weakly 
bound to the island because they have a smaller number of bonds, so minimising their number 
is thermodynamically favourable. This leads to large and circular islands being the most stable 
configuration for the system because it minimises the interface free energy. 



In Fig. 5.1, larger temperatures present this thermodynamically favourable state. This is 
caused by the high thermal energies that lead to large adatom mobility and the system can easily 
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Figure 5.1: Physical lattice at temperatures (a) T = 550 K, (b) T = 600 K, (c) T = 650 K, (d) T = 700 K. In 
all cases, E a = 1.00 eV, Et, = 0.20 eV, F = 1 ML/s and 6 = 0.1 ML. Black squares represent occupied sites and 
white squares unoccupied sites. 



reach the most favourable configuration. Under these conditions, islands that are not very large 
can break up and only a few nucleations will eventually lead to stable islands, which will be 
fewer in number. Also, the islands will be larger in size because there will be less competition 
for adatom incorporation due to the presence of a smaller number of competing islands. 



When the flux is turned off all the different cases depicted in Fig. 5.1 will tend to configura- 
tions minimising the free energy for large times, even when the thermal energies do not facilitate 
it. A common mechanism that leads to coarsening in submonolayer epitaxial systems is Ostwald 
ripening [3lJ [32] . 
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Figure 5.2: Adatom density n and island density N as a function of time. The system corresponds to the 
standard KMC model with bond energy E b = 0.20 eV and T = 650 K. 



These features can be studied quantitatively by looking at the island size distribution and 
at the time and temperature dependence of adatom and island densities as described in the 



following subsections. For clarity, we plot in Fig. 5.2 the time dependence of both adatom 
and island densities for the parameters Ef, = 0.20 eV and T = 650 K. Such profiles are used 
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in the various analyses presented below for the island size distribution, the adatom density at 
nucleation n nuc , the adatom density at equilibrium n eq and the island density at equilibrium 



N eq . In the instances shown in Fig. 5.2 the deposition flux is turned off at a coverage of 9 



0.10 ML (corresponding to t = 0.10 s for F = 1 ML/s) and the system is then allowed to relax 
for a further 0.05 seconds. The main features discussed in the previous chapter for such profiles 
are clearly observed. 

5.2.2 Island size distribution 

The island size distribution is a quantity readily available in KMC studies, but difficult to obtain 
in mean field approaches such as rate equations. 
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Figure 5.3: Island size distribution for the standard KMC model with bond energy E — 0.20 eV. 



In Fig. |5.3| we plot the island size distribution as it is at a coverage of 9 
energy E h 



0.1 ML, for bond 



0.20 eV and for a range of temperatures. The distributions are calculated when 
the system has been allowed to relax after turning the deposition flux off. The island sizes are 
normalised in such a manner that the area under the curve is unity, and the horizontal axis 
represents s/s av for island size s and average island size 



EsNs 



(5.2) 



where N q is the number of islands of size s. 



It can be seen in Fig. 5.3 that as temperature increases the distribution becomes narrower. 
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At large temperatures, islands can break up more easily than at low temperatures, so we expect 
that there is some correction in island sizes and all islands approach an ideal size. This analysis 



quantifies the discussion above based on snapshots of the lattice in Fig. 5T 

To exemplify this island size correction at high temperatures consider a system where two 
islands form nearby. At low temperatures, the islands will grow by attachment of adatoms, 
but because they are nearby they compete for free adatoms in the surface, so both islands 
will in general be smaller than the average island size. However, the same situation at high 
temperatures can result in the break-up of one of the two islands, and the incorporation of its 
adatoms into the other, so that only one island is left, and its size is closer an average size. 
Note that in the same manner, if an island is formed in a region where it is isolated, for low 
temperatures it will incorporate all adatoms in its vicinity and thus grow to larger than average 
sizes. However, at large temperatures the detachment of adatoms from the island will be more 
significant and the size of the island will again be closer to that of others. 

The data for island size distribution is not as clear as the data for adatom or island concen- 
trations, specially for large temperatures. This is because for high temperatures very few islands 



are formed on the lattice, in Fig. 5.1 only 10 islands for T = 700K. Furthermore, the islands 
formed are larger in size, so the range of possible island sizes is broader. These two characteris- 
tics mean that to obtain useful data for island size distributions we need to calculate an average 
over many different instances of the system evolution. For instance, in the case depicted in Fig. 



5.3 for T = 700K, up to 10,000 instances were used. 



5.2.3 Adatom density 

The experimental results reported by Loginova and co-workers O 13] represent some of the 
cleanest data on adatom density as a function of time for epitaxial growth experiments. As 
described above, the adatom density presents interesting features, and we will concentrate on 
its values at the onset of nucleation, n nuc , and at equilibrium with islands, n eq . 

The behaviour of n nuc is highly dependent on the strength of the bond between adatoms 



as shown in Fig. 5.4 For low bond energy the critical island size is large, that is, the size 
required for a cluster of atoms to eventually become an island is large because small clusters will 
break-up with large probability due to the weakness of the bonds keeping the atoms together. 
Furthermore, as temperature increases bonds break more often, so we expect that with increasing 
temperature nucleation occurs at larger adatom concentrations. In contrast, for large bond 
energies the thermal vibrations cannot easily break the bonds, and the critical island size is 
reduced to a few atoms even for high temperatures. Then the limiting process for nucleation 
is adatom diffusion because as soon as adatoms become nearest neighbours nucleations occur. 
Therefore, as temperature is increased, the larger diffusion of adatoms means that the adatom 
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Figure 5.4: Adatom nucleation density n nuc as a function of temperature for nearest neighbour bond strengths 
E b = 0.10, 0.20, 0.30 eV. 



nucleation occurs sooner and n nuc decreases with temperature. 



In Fig. 5.4 we plot n nuc as a function of temperature, and for different adatom bond energies 
Ef,. Each data point has a statistical error associated with it that is of order ~0.01 of the value 
plotted so it is not depicted. We can observe both the low bond energy limit with E\, = 0.10 
eV and the large bond energy limit with E\, = 0.30 eV. With bond energy E^ = 0.20 eV we 
can see an intermediate scenario where the lower temperature range is diffusion dominated and 
the adatom nucleation concentration decreases with temperature, but at larger temperatures 
island break-up becomes dominant and a larger adatom concentration is required for nucleation 
because the critical island size increases. Note also that as expected, for lower bond energies 
n nuc is much larger than for large bond energies. 



The adatom density at equilibrium is shown in Fig. 5.5 , plotted as a function of temperature 
and for different E^. In all cases the equilibrium concentration increases with temperature 
because detachment from islands dominates at larger temperatures. For low binding energy 
Eb = 0.10 eV, n eq is large because detachment rates are very large. In contrast, for higher bond 
energies the adatoms do not detach from islands with such large rates and n eq almost vanishes. 
For the lower temperatures and higher binding energies equilibrium is only reached after a long 
time because low diffusion over the substrate means adatoms do not attach to islands frequently, 
but strong bonds lead to large islands, so the system takes a long time to equilibrate. 
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Figure 5.5: Adatom equilibrium density n eq as a function of temperature for nearest neighbour bond strengths 
E b = 0.10, 0.20, 0.30 eV. 



5.2.4 Island density 

In this subsection we study the equilibrium island density N eq for the standard model. Exper- 
imental data on island density is available for the graphene system, and the RE approach fails 
to describe it. Therefore, a careful analysis of island density within KMC could deepen our 
understanding of the growth kinetics of graphene. 

In the bond-counting scheme used for the standard model, a critical island size is not well- 
defined. Therefore, the definition of island is not very clear. For instance, for low bond energies 
between nearest neighbours, dimers are very unstable and break-up frequently. However, for 
large bond energies, dimers are very stable and lead to large islands in most cases. In the 
standard model, clusters of size larger than 10 are treated as islands. This choice might not be 
strictly accurate, but the trends to be observed in the data arise clearly, therefore simplicity has 
prevailed in the choice. 



In Fig. 5.6 we plot the equilibrium island density as a function of temperature. For the 
stronger bond cases Ef, = 0.20, 0.30 eV, a clear decreasing trend for N eq with temperature is 
observed. This is in agreement with the physical lattices above. Note also that the concentration 
of islands is larger for the case with stronger binding. This is caused by the higher difficulty of 
small clusters to break up if bonds are strong, and therefore more nucleations are expected in 
this regime. The low bond energy case E^ = 0.10 eV presents a different behaviour. For the 
temperature increase from 550 K to 600 K, a decrease in N eq is observed as expected. However, 
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Figure 5.6: Equilibrium island density N eq as a function of temperature for nearest neighbour bond strengths 
E b = 0.10, 0.20, 0.30 eV. 



further temperature increases do not lead to a clear decrease, instead the island concentration 
remains approximately constant. This behaviour is very particular of the low bond energy 
case where islands cannot grow very large due to the high adatom detachment rates. The 
observed behaviour of island density means that already at T ~ 600 K the adatom density on the 
substrate determines the minimum island density attainable with the high adatom detachment 
rate present. Then, further temperature increases do not lead to a smaller island concentration. 
This behaviour is confirmed by the physical lattice in this regime shown in Fig. |5.7[ where 
temperatures higher than T = 600 K do not change the aspect of the system as significantly as 



for the case depicted above in Fig. 5.1 The islands decrease slightly in size due to the higher 
temperatures causing more frequent bond breaking. 

The behaviour of the system suggests that we use a nearest neighbour bond energy of the 
order Eb ~ 0.20 eV or above for further studies of graphene. This is because lower binding 
energies will not lead to growth of large islands. 



5.2.5 Error analysis 

The different values presented above have a statistical uncertainty associated with them due to 
the stochastic nature of the KMC method. These statistical errors are minimised by averaging 
the results over a large number of experiments, and this has been done in a manner that the 
errors of a given quantity are about 1% of the value of the quantity throughout this work. 



26 




Figure 5.7: Physical lattice at temperatures (a) T = 550 K, (b) T = 600 K, (c) T = 650 K, (d) T — 700 K for 
a bond energy of E b — 0.10 eV. Black squares represent occupied sites and white squares unoccupied sites. 



These small statistical errors are not shown in the plots because they are too small to be seen 
appropriately. In the quantity N eq the statistical errors can be up to 10% of the value quoted, 
but even in these circumstances the trends discussed can be clearly observed. 

Another error consideration comes in when evaluating the correctness of the model and its 
implementation. The model described here is a simple bond-counting epitaxial growth model 
that is widely used in such studies. Its implementation has been partially checked by finding 
the expected results for such bond-counting model. In the next chapter we consider this point 
further for the KMC model to describe epitaxial graphene growth. 
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Chapter 6 



Epitaxial graphene growth kinetic 
Monte Carlo model 

In this chapter we present a KMC model for the study of the kinetics of epitaxial graphene 
growth. The experimental observations reported about the system [31 H] are explained with the 
model, and the relevant physical processes determining the observed behaviour are identified. 

However, direct comparison with experiment is not possible for two reasons. First, the values 
of the energy barriers associated with the different processes are not available yet. The present 
work determines the relevant physical processes to include in the study of epitaxial graphene 
growth. Therefore, the next step is to determine the exact energetics of these identified processes, 
for instance using ab initio techniques. Second, the deposition flux used experimentally is very 
low, resulting in time scales that are computationally inaccessible to the present KMC model. 

6.1 Description 

We introduce a KMC model for epitaxial graphene growth in a 200x200 square lattice with 
periodic boundary conditions. We include adatoms, islands and an intermediate species formed 
by four adatoms arranged in a square which we call tetramers. Experiments indicate that 
the intermediate species in the graphene system are 5-atom clusters, but for convenience with 
the underlaying square lattice we choose the simpler symmetric tetramer. The presence of 5- 
clusters in graphene growth is most probably determined by the stability of such structures on 
the substrates (observed in grand canonical Monte Carlo simulations for Niquel Ni(lll) [33]), 
but we postulate that the kinetics are determined only by the presence of this intermediate 
species rather than by its details. This justifies the use of tetramers which represent the stable 
symmetric structure in our model. The validity of this assumption can be confirmed later with 
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the results obtained. In a similar manner, the use of a square lattice rather than a triangular 
lattice can be justified for simplicity of the model. 

Deposition occurs as described above for the standard model via individual atoms. Exper- 
imentally, both carbon atoms and ethylene molecules are deposited, but the behaviour of the 
growth system is independent of the deposited species [I] and for simplicity we only consider 
individual atom deposition. 

Adatoms diffuse over the substrate in the same manner as for the standard model. They only 
interact to form tetramers and to attach to islands. A tetramer is formed whenever four adatoms 
come together in a square configuration on the lattice. From that moment onwards, the tetramer 
is treated as an individual species formed by four occupied sites in a square configuration. 
Attachment of adatoms to islands occurs as described above for the standard model whenever 
an adatom diffuses to a site with nearest neighbours belonging to an island. After an adatom has 
attached to an island, it is still treated as an individual atom, but subject to the bond-counting 
scheme described above. We stress this only occurs for adatoms that diffuse to sites with nearest 
neighbours belonging to islands: if a tetramer is a nearest neighbour, there is no interaction with 
the adatom based on the fact that tetramers are a very stable and independent species. 




tetramer 



Figure 6.1: Schematic of possible nucleation configurations for (a) j — 2, (b) j = 4 and (c) j = 6, where j is 
the number of tetramers. For the asymmetric cases j = 2, 6; a n/2 rotation with respect to the cases depicted is 
also a valid nucleation configuration. Note that in the diagram a black circle represents four adatoms in a square 
configuration and a square lattice site represents four square lattice sites in the actual simulation. 



Tetramers can diffuse over the substrate, in a manner analogous to adatom diffusion. The 
only difference is that the four adatoms forming the tetramer move at the same time and in the 
same direction as an individual unit. Tetramers only interact with other tetramers to nucleate 
and with islands to attach to them. Nucleation occurs when a number j of tetramers come 



together on the lattice in a symmetric manner as indicated in Fig. 6.1, where three possible 
values for j are shown. As soon as a nucleation occurs, the tetramers lose their identity and 
the adatoms forming them become members of an island and are subject to nearest neighbour 
interactions via the bond-counting scheme. Tetramers can also attach to existing islands, and 
such events occur whenever the two atoms of a side of the tetramer become nearest neighbours 
to sites occupied by islands, in such a manner that a side of the tetramer square is in contact 
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with an island. This is schematised in Fig. 6.2 When a tetramer attaches to an island, it again 
loses its identity and the individual atoms become subject to nearest neighbour interactions. 



adatom belonging to tetramer 
adatom belonging to island 



Figure 6.2: Schematic of tetramer diffusion and attachment to islands. 



We tune the temperature whilst keeping the energy barriers shown in Table 6.1 at fixed 
values. Note the tetramer diffusion barrier is kept lower than the adatom diffusion barrier to 
promote tetramers to the dominant species in graphene growth kinetics. 



Description 


Parameter 


Value 


Coverage 


9 


0.10 ML 


Adatom diffusion barrier 


E a 


1.00 eV 


Tetramer diffusion barrier 


E t 


0.80 eV 


Adatom bond energy 


E b 


0.20 eV 



Table 6.1: Values of the basic parameters used in the KMC model to describe epitaxial graphene growth 



The model presented here belongs to the generic models class. A possible variation of this 
model for the simulation of epitaxial graphene growth is described in Appendix [Aj 

In the following section we investigate how different processes superimposed on top of this 
basic model determine the rich behaviour observed experimentally for the kinetics of epitaxial 
graphene growth. 



6.2 Results and discussion 

In this section we study the three experimental observations by Loginova and co-workers on the 
epitaxial graphene growth system [31 E]- For convenience, we repeat them here: 

1. The adatom density at the onset of nucleation increases with increasing temperature. 

2. The adatom density at equibrium n eq is roughly half of the adatom density at the onset 
of nucleation n nuc , namely, n nuc ~ 2n eq . 

3. The island density at equilibrium decreases rapidly with increasing temperature. 
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In the following subsections we study them individually. There is only preliminary data for the 
island size distribution, and it is presented in Appendix [B| 



6.2.1 Adatom density at nucleation 



It was discussed in Section |3.1| that the adatom density at nucleation, n nuc , is experimentally 
found to increase as a function of temperature for the graphene growth system. In standard 
growth systems, we found in Chapter [5] that n nuc could present different behaviours with tem- 
perature depending on the binding energy between adatoms. In this section we will investigate 
this effect on the KMC graphene model and find a more complex behaviour. 

Inspired by experiments (23J El] and theoretical studies [5], it appears that for graphene 
the critical nucleus sizes are of the order of tens of adatoms. In our KMC model we have 
incorporated this observation by allowing nucleations to occur only when a certain number j of 
tetramers come together in specific configurations. Because islands only start appearing with 
these large sizes, then the analysis on the standard model with low binding energy might not 
be applicable for graphene. This is because even for low bond energies, the incipient islands 
are so large that dissociation might be very unlikely and the growth process would be diffusion 
limited. Such a scenario would lead to a decrease in n nuc with increasing temperature, even for 
low energy bonds. Furthermore, it was concluded in the previous chapter that binding energies 
which are too low do not result in the appearence of large islands. 

Even though we might have a diffusion limited growth system, epitaxial graphene growth is 
dominated by an intermediate stable species, and this difference with other growth systems can 
explain the increase of n nuc with temperature, even with the presence of relatively large incipient 
islands. Physically, a process is needed by which the concentration of clusters giving nucleation 
is decreased with increasing temperature to compensate for the larger cluster formation and 
diffusion rates. Such a process could be temperature dependent cluster break-up, included in 
the KMC model with a rate in the usual Arrhenius form subject to a cluster break-up energy 
barrier E^. 



In Fig. 6.3 we show the temperature dependence of n nuc for different nucleation sizes j and 
for a range of cluster break-up energies Ej.. We also plot the corresponding quantities for the 
standard model discussed above. 

For j = 2 we find that for all cluster break-up energy values, n nuc decreases with increasing 
temperature. This indicates that, as discussed above, an incipient island formed by N = 8 
adatoms is large enough not to dissociate, so that a low binding energy between adatoms cannot 
lead to an increase of n nuc with temperature as for the standard model even at high temperatures. 
Furthermore, cluster break-up is also found not to be sufficient to appropriately describe the 
temperature behaviour of n nuc of graphene growth. Even though the tetramer population is 
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Figure 6.3: Adatom density at the onset of nucleation as a function of temperature for (a) standard model, 
(b) nucleation with j = 2, (c) nucleation with j = 4, (d) nucleation with j = 6. In the standard model different 
nearest neighbour bond energies are considered, and in the tetramer model different tetramer break-up energies 
are considered. 



reduced as temperature increases, for j = 2 nucleation is still a frequent event and the break-up 
rate is not high enough to act as a significant deterrent for nucleation. This means the system 
behaves as a typical diffusion limited growth system and n nuc decreases with temperature. 

For j = 4 more interesting behaviour starts appearing. For zero or very small cluster break- 
up, we still find a decreasing trend in n nuc . Again, this is due to the behaviour as a typical 
diffusion limited growth system. However, as is decreased, an increase in n nuc is observed 
for the larger temperature ranges. With j = 4, nucleation events are less likely than with j = 
2, so at large temperatures where cluster break-up is a significant effect, it can counteract the 
larger tetramer formation and diffusion rates. For a critical concentration of tetramers needed 
for nucleation, the high tetramer break-up rate means that larger adatom concentrations are 
needed to lead to the necessary tetramer concentration for nucleation. This leads to the increase 
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of 

n nuc with temperature. It is not observed for low temperatures because diffusion still prevails 
over tetramer break-up in that regime. 

For j = 6, with the appropriate choice of Ek, an increase of n nuc with temperature is 
observed in the entire temperature range. This is the experimentally observed behaviour of 
epitaxial graphene growth. The large value of j means that the tetramer concentrations needed 
for nucleation are so large that tetramer break-up can play a significant role even for low tem- 
peratures. From now on and unless otherwise stated we are going to use j = 6 and E^ = 1.40 
eV for all further studies. 

The above reasoning depends on the fact that the density of clusters at the onset of nucleation 
is independent of the cluster break-up energy barrier. This is indeed the case as shown in Fig. 



6.4 for j = 4, 6 and similar behaviour occurs for j = 2. 
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Figure 6.4: Tetramer density at the onset of nucleation as a function of temperature for (a) nucleation with j — 
4 and (b) nucleation with j = 6. 

In agreement with the rate equations study by Zangwill and Vvedenksy [5j, both cluster 
break-up and large incipient island sizes are found to be necessary to explain the increase of 
n-nuc with temperature. 



6.2.2 Adatom density at equilibrium 

Adatom density at equilibrium is experimentally found to be roughly half of the adatom density 
at the onset of nucleation, n nuc ~ 2n eq . Such large n eg require the presence of a large attachment 
barrier for adatoms to graphene. The presence of such a barrier initially prompted the proposal of 
a 5-cluster attachment mechanism by Loginova and co-workers [3] . In the KMC model presented 
in this thesis, two physical processes determine the predominant growth of graphene by cluster 
attachment. The first one is a larger diffusion for tetramers than for adatoms, and the second is 
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the presence of an adatom attachment barrier. It is only the second that determines the adatom 
concentration at equilibrium, and this effect is discussed here. 

We consider a temperature dependent attachment barrier for adatoms to graphene. Com- 
putationally, it is incorporated in the KMC model via a Boltzmann factor e~ Ep / kBT such that 
when an adatom moves to a site with nearest neighbours belonging to an island, it will attach 
to the island with probability p < e ~ E vl k BT _ This condition makes adatom attachment more 
difficult at low temperatures than at high temperatures. Tuning the energy barrier E p results 
in different values for n eq . 
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Figure 6.5: Ratio of equilibrium adatom concentration n eq over adatom concentration at the onset of nucleation 
Tinuc as a function of temperature. Different adatom attachment barriers are considered. 



In Fig. 



6.5 



we plot the ratio n eq /n nuc as a function of temperature for a range of adatom 
attachment barriers. We calculate n eq from the KMC model by turning the deposition flux off 
and allowing the system to relax, that is, by allowing the system to reach equilibrium between 



adatoms, tetramers and islands. As can be seen in Fig. 6.5, the value of E p tunes the adatom 

can be reached in our 



In 



eq 



equilibrium concentration, and typical experimental values n nv 
model for E p = 0.30 eV. The ratio is also seen to be roughly constant for the temperature range 
studied, although there is a slight increase with temperature. From now on unless otherwise 
stated we are going to take E p = 0.30 eV. 

Recalling n eq for the standard model in the previous chapter, we find significant differences. 
For the case with adatom bond energy of Ej, = 0.20 eV as used here, the standard model 
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resulted in vanishingly small equilibrium concentrations. Therefore, it is clear that an adatom 
attachment barrier is essential in obtaining the correct behaviour. 

The RE approach by Zangwill and Vvedensky |5j also describes n eq accurately. However, in 
their model there is no adatom attachment barrier, so the only process that can also explain 
this behaviour is a large adatom detachment rate from islands. The lack of spatial information 
in the RE means that the capture zone of islands is smaller than in the KMC model, so adatoms 
attach to islands less frequently. This leads to higher values for n eq in the RE than in KMC 
models without an adatom attachment barrier. 

6.2.3 Island density at equilibrium 

Island density at equilibrium is experimentally found to decrease rapidly with increasing tem- 
perature. In the KMC model used above for the description of graphene growth, as soon as a 
nucleation occurs, tetramers are more likely to attach to the existing island than to give rise 
to further nucleations. Therefore, to have a significant number of nucleations, the density of 
tetramers needs to be kept high by some mechanism. In order to try to describe this, we incor- 
porate a tetramer attachment barrier E q in the same manner as an adatom attachment barrier 
above. 
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Figure 6.6: Equilibrium island concentration N eq as a function of temperature. Different tetramer attachment 
barriers are considered. 
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In Fig. 6.6 we plot N eq as a function of temperature for a range of tetramer attachment 
barriers E q . It can be seen that, by tuning the tetramer attachment barrier, a large decrease in 
island density can be engineered into the model. For the standard model discussed above, N eq 
decreased slowly as a function of temperature. For instance, for a bond energy E\, = 0.20 eV, 
an increase of 100 K gave a decrease by a factor of 4. In the present case, an increase of 100 K 
for E q = 0.30 eV gives an order of magnitude decrease in N eq . 

Note that in order to get a large decrease in island density at nucleation, tetramer attachment 
barriers of the same order as adatom attachment barriers need to be used. In this scenario, the 
only property of the system that makes tetramer attachment a predominant effect in graphene 
growth is the larger diffusion of tetramers as compared to adatoms. In the next section we 
are going to confirm that under these circumstances tetramer attachment is still the dominant 
effect. 



6.3 Error analysis: island growth velocity 

In this section we present an error analysis on the KMC model described in this thesis. Statistical 
errors due to the stochastic nature of the method are already discussed in the context of the 
standard model above, and that discussion is not repeated here. Instead, we concentrate in the 
evaluation of the model used and its computational implementation. We look at the velocity of 
island growth and find that the KMC model is consistent with the m-cluster model proposed 
experimentally for the growth of epitaxial graphene. 

In Ref.[3j the island growth velocity is found to be a nonlinear function of the adatom 
concentration, ultimately leading to the proposal of a growth model for graphene with a 5- 
cluster intermediate species. As discussed in Section [3. 1[ postulating that the energy barrier for 
adatom attachment to islands is larger than the energy barriers for clusters of size m to form 
and later to join islands, it can be shown that the island growth velocity obeys 



D 



n , 

1 



(6.1) 



where n is the adatom concentration, n eq is the adatom concentration at equilibrium and B is a 
temperature dependent constant. Experimentally, it is found that m ~ 5, and the intermediate 
species are clusters formed by 5 carbon atoms. In our KMC model, where tetramers play the 



role of the intermediate species, we expect to have a power law with exponent m ~ 4 in Eq.(6.1 ) 
This result will provide a consistency check for the model used in this work. 

The island growth velocity v is defined as [3] 
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for island perimeter P and area A. For simplicity, we consider circular islands with P = 2ttR 
and A = irR 2 for radius R. The circular island approximation is expected to be accurate at 
high temperatures where we have seen that the system tends to large circular islands trying to 
minimise the edge free energy. With a total area A = N for iV carbon atoms in the island and 
an island radius R = \J N/tt, then 

P = 2v / vriV and A = N. (6.3) 

The estimate for the island growth velocity v(t) at time t in the circular island approximation 
is then 

. . 1 AA 1 N(t + At)-N(t) . 

v(t)~—— ~ — — ^ — . (6.4) 

W P At 2yGfN{t) At v ' 
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Figure 6.7: Island growth velocity as a function of adatom coverage at T — 700K for tetramer attachment 



barriers (a) E q = 0.10 eV and (b) E q = 0.30 eV. A fit to Eq.(6.1 \ is also shown 



In Figure 6.7 we plot the v as a function of n for tetramer attachment barriers E q = 0.10, 0.30 



eV. The other parameters are those in Table 6.1 and = 1.40 eV, E p = 0.30 eV. The data 
corresponds to T = 700 K, the large temperature limit discussed above, and it is then found 
that n eq = 0.032 ML for E q = 0.10 eV and n eq = 0.039 ML for E q = 0.30 eV. 



The fit has been done using the Mathematica package with the model given by Eq.(6.1 ) and 



leaving as free parameters to fit B and m. The least-squares fit parameters are shown in Table 



6.2 where the quoted errors are the standard error. We note that the fit parameter m is highly 
dependent on the value of n eq which needs to be determined with high accuracy. 

In both cases the fit value for m indicates that the implementation of the model is consis- 



ts 



E„ = 0.10 eV E„ = 0.30 eV 



B 299.6 ± 88.3 
m 4.09 ± 0.43 



94.9 ±43.7 
4.44 ±0.80 



Table 6.2: Least-squares fit values of the parameters B and m in Eq.(6.1l for tetramer attachment barriers 
E a = 0.10 eV and E a = 0.30 eV. 



tent with the m-cluster model proposed by Loginova and co-workers, in our case with m ~ 4 



because we used tetramers as the intermediate species. In the model described by Eq.(6.1), the 



attachment to graphene is via m-clusters. In the present study, the growth is via both tetramers 
and adatoms so we expect some discrepancies with Eq.(6.1). For E q = 0.10 eV, m-cluster at- 



tachment is more predominant over adatoms than for the case with E q = 0.30 eV, and indeed 



we observe a better fit. However, the discussion above in Subsection 6.2.3 about N eq suggests 
that the appropriate value for tetramer attachment barrier is E q = 0.30 eV. Even though the fit 
is not as accurate as for E = 0.10 eV, the velocity analysis here suggests that E q = 0.30 eV is 



still consistent with Eq.(6.1). This indicates that the larger tetramer diffusivity still makes this 



species the predominant one in graphene growth. 

Computationally, island sizes at time intervals of At = 2.5 ms have been used to obtain the 



results shown in Fig. 6.7 Large islands can break-up or coalesce during growth, and then v is 
not evaluated as the values do not correspond to growth by attachment of adatoms or tetramers 
only. Also, a deposition flux F = 10 ML/s has been used to explore the larger values of the 
density n. 

Different approaches to check the validity of the model and its implementation can be con- 
sidered outside of the results obtained within the model. In Appendix [C] we present a RE 
comparison with the present KMC model that further validates it. 
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Chapter 7 

Conclusions 



7.1 Summary 

In this thesis we have introduced a KMC model to study the physical processes governing the 
kinetics of epitaxial graphene growth on metal substrates. All the experimental observations 
have been explained within the model, and the relevant kinetic paths determining them have 
been identified. 

The behaviour of adatom density at the onset of nucleation has been explained by the 
inclusion of both a large critical island size and a tetramer break-up rate. These conclusions are 
in agreement with previous RE studies. 

The adatom density at equilibrium has been found to be highly dependent on the adatom 
attachment barrier to islands. With no such barrier, the adatom density at equilibrium is van- 
ishingly small, but tuning the energetics of such a process can lead to the desired concentration 
of adatoms at equilibrium. 

The island density at equilibrium has been found to be dependent on the attachment barrier 
of tetramers to islands. Without such barrier, few nucleations occur because tetramers are more 
likely to attach to an existing island than to nucleate new islands. Only with the inclusion of a 
large barrier the desired behaviour in the island density at equilibrium is obtained. This means 
that the dominance of tetramers in graphene growth is via their larger diffusivity rather than 
due to larger attachment barriers for adatoms than for tetramers. 

7.2 Further work 

The present work has allowed the determination of the relevant physical processes governing the 
particular behaviour of the kinetics of epitaxial graphene growth. The next natural step is to use 
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this knowledge to construct an accurate KMC model to be compared directly with experiment. 
First, the various energy barriers in the problem need to be determined. A standard way of 
accomplishing this is by electronic structure calculations using well-established first-principles 
methods such as density functional theory. Also, both the metal substrate (lattice) and the 
graphene structure need to be hexagonal. Finally, an exact intermediate 5-cluster species instead 
of the more convenient tetramers needs to be included. Difficulties in this model might arise 
when comparing with experiment because high temperatures and large time intervals need to be 
studied, and this might be computationally too demanding. The code can be further optimised, 
but some approximations might still be necessary to reach the scales of interest. 

The detemination of the relevant energy barriers for epitaxial graphene growth by means of 
ab initio techniques can also be used in similar studies such as lattice- free KMC. 



Words in text: 9986. 
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Appendix A 

Island growth model 



In this appendix we consider the role of islands in the KMC model presented in Chapter [6j 

The model presented in this thesis falls within the generic model class for the treatment of 
islands, but the rules on the existence of tetramers and nucleation of islands make it different 
from a pure bond-counting model. Generic models [30] are those with local rules such as bond 
counting according to which the system evolves, and are different from the so-called tailored 
models. In the latter, some specifications on the system evolution are made based on the 
knowledge already existing about the system. For instance, a critical island size can be defined 
above which attachment of adatoms to islands is irreversible. Another possibility is the prescribed 
island growth sequence model [30] that we will discuss here in the context of graphene. 

In a prescribed island growth sequence model, adatoms that attach to an island are instan- 
tenously moved to a prescribed position on the island independently of the attachment position. 
This is equivalent to having an infinite edge diffusitivity to the most stable sites. It is used in 
systems known to grow in certain shapes, for instance in growth on metal (100) where the islands 
are known to be nearly squares and each new attached adatom is immediately incorporated in 
the sequence to give the desired square shape. 

For graphene growth via 5-clusters, the structure of the hexagonal lattice makes the 5- 
clusters have the appropriate size to represent an attachment in which growth only occurs via the 
formation of new complete hexagons. In Fig. |A.1| we can see such configurations schematically, 
where we can note that 3-clusters would also lead to such hexagon by hexagon growth and 
5-clusters actually lead to double hexagon growth. This indicates that graphene grows in a very 
prescribed manner, and the tetramer attachment and bond counting scheme used in this thesis 
might not be entirely accurate in describing it. A prescribed island growth sequence, by which 
tetramer attachment leads to an immediate relaxation into a more stable configuration might 
be more appropriate. 
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The island shape observed in graphene and shown in Fig. A. 2 also suggests a very specific 
growth sequence, possibly the result of this attachment mechanism. 
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Appendix B 

Island size distribution 



In this appendix we present the preliminary data for the island size distribution of the KMC 
model for the description of epitaxial growth. For the reasons discussed in Chapter [5j data 
for island size distribution is very noisy, and to get significant results many instances of the 
system need to be used. Time constraints have limited the results to the case of T = 500 K for 



the graphene model so far. These are shown in Fig. B.l for the parameters in Table 6.1 and 
Ek = 1.40 eV, E p = 0.30 eV and E q = 0.30 eV as usual. The shape observed is to be expected 
for the island size distribution. However, any detailed analysis will require the behaviour of 
island size distribution with temperature. 
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Figure B.l: Island size distribution for the graphene KMC model at T = 500 K. 
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Appendix C 

Rate equations and kinetic Monte 
Carlo comparison 



In this appendix we compare the KMC model presented in this thesis with the RE approach by 
Zangwill and Vvedensky j5]. The objective of the comparison is to validate the model used in 
this work. 





RE 


KMC 


E a 


1.00 eV 


1.00 eV 


E t 


0.80 eV 


0.80 eV 


Ek 


1.40 eV 


1.40 eV 


E b 




0.20 eV 


Ei 


1.60 eV 





Table C.l: Least-squares fit values of the parameters B and m in Eq.(6.1l for tetramer attachment barriers 



E q = 0.10 eV and E q = 0.30 eV. 



The rate equations model used is the one presented by Zangwill and Vvedensky [5] as de- 



scribed in Sectin 3.2 above. This RE model includes parameters to describe adatom and tetramer 
diffusion, tetramer break-up and adatom detachment E{ from islands. The KMC model has cor- 
responding parameters for the first three quantities, but for the latter it is the adatom bond 
energy between adatoms in an island that maps the effect. In order to carry out the comparison, 
RE have been used with i = 4 and j = 6 and KMC has been used with E p = E q = 0. The rest 



of parameters are shown in Table C.l A temperature of T = 650 K is used with a flux F = 1 
ML/s. The total coverage reached is = 0.2 ML corresponding to 0.2 seconds, and the system 
is allowed to relax for a further 0.3 seconds. 



Figure C.l shows the adatom concentration n and the tetramer concentration c as a function 



of temperature. In both cases, KMC and RE coincide up to the nucleation point at about 
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(a) Adatom density (b) Tetramer density 

Figure C.l: Comparison between RE and KMC model for (a) adatom density and (b) tetramer density. 



t ~ 0.1 s, but a clear disagreement appears after this. The RE overestimate the concentration 
of both adatoms and tetramers. In the KMC approach, islands have spatial extent, so their 
adatom and tetramer capture zones are large. However, in the RE the islands are point-like, 
and their capture zones are reduced, thus leading to a smaller decrease in the concentration of 
adatoms and tetramers. 
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Figure C.2: Comparison between RE and KMC model for the island density. 



In Fig. C.2 we show the comparison between RE and KMC for the island density. Unlike the 



above, this quantity is badly reproduced by the rate equations as already discussed in the work 
by Zangwill and Vvedensky. The island concetration is underestimated by the rate equations, 
and again this is attributed due to the lack of spatial information [5] . 
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